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Abstract 



We have observed 7 new transits of the 'hot Jupiter' WASP-5b using a 61 cm 
telescope located in New Zealand, in order to search for transit timing variations 
(TTVs) which can be induced by additional bodies existing in the system. When 
combined with other available photometric and radial velocity (RV) data, we find 
that its transit timings do not match a linear ephemeris; the best fit x 2 values is 32.2 
with 9 degrees of freedom which corresponds to a confidence level of 99.982 % or 3.7 o. 
This result indicates that excess variations of transit timings has been observed, due 
either to unknown systematic effects or possibly to real TTVs. The TTV amplitude 
is as large as 50 s, and if this is real, it cannot be explained by other effects than 
that due to an additional body or bodies. From the RV data, we put an upper limit 
on the RV amplitude caused by the possible secondary body (planet) as 21 m s" 1 , 
which corresponds to its mass of 22-70 M e over the orbital period ratio of the two 
planets from 0.2 to 5.0. From the TTVs data, using the numerical simulations, we 
place more stringent limits down to 2 M® near 1:2 and 2:1 mean motion resonances 
(MMRs) with WASP-5b at the 3 a level, assuming that the two planets are co-planer. 
We also put an upper limit on excess of Trojan mass as 43 M ffi (3 a) using both RV 
and photometric data. We also find that if the possible secondary planet has non- or 
a small eccentricity, its orbit would likely be near low-order MMRs. Further follow-up 
photometric and spectroscopic observations will be required to confirm the reality of 
the TTV signal, and results such as these will provide important information for the 
migration mechanisms of planetary systems. 

Key words: stars: planetary systems: individual (WASP-5) — techniques: pho- 
tometric 

1. Introduction 

More than 400 extrasolar planetary systems have been found to date by several tech- 
niques, including pulsar timing, radial velocity (RV), transit, microlensing, direct imaging, and 
astrometry. Among them, more than 40 systems have been revealed to have multiple planets, 
most of which have been detected by the RV technique (Wright et al. 2009). Considering obser- 
vational limits, multi-planetary systems are thought to be common — that is naturally expected 
from the standard planetary formation mechanism, called core accretion models (e.g., Hayashi 
1981; Pollack et al. 1996; Kokubo & Ida 2002). Increasing the number of multi-planetary 
systems and studying their nature are important to improve our understandings of planetary 
formation mechanisms and dynamics of planetary systems. 

Searching for transit timing variations (TTVs) from a constant period ephemeris is 
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another method for identifying multi-planetary systems. If a transiting planet is the one and 
only body orbiting its host star, its orbital period should be constant. On the other hand, if 
another perturbing body exists in the system, the orbital period will no longer be constant. 
Therefore we can find additional planets by probing TTVs, even if they do not transit their 
host star (Agol et al. 2005; Holman & Murray 2005). Furthermore, if they also transit, one can 
solve for their mass and orbital elements without RV measurements — at least in principle. 
Not only additional planets but also exomoons (Kipping 2009) or Trojans (Ford & Holman 
2007) can be detected via the TTV method. TTVs can also be caused by other effects, such 
as orbital precession due to effects of the general relativity, spin-induced stellar quadrupole 
moment, and tidal deformations of the star and the planet (e.g., Iorio 2010, and references 
therein), as well as variations in the quadrupole moment driven by stellar activity known as 
the Applegate effect (Applegate 1992; Watson & Marsh 2010) and orbital decay due to tidal 
dissipation (Hellier et al. 2009). However, the impact of these effects on the orbital period is 
relatively tiny compared to that due to additional orbiting bodies. 

The TTV method is especially sensitive to planets in orbits near mean motion reso- 
nances (MMRs) with transiting planets. Even a sub-Earth mass planet in low-order MMR is 
detectable using relatively small-size (0.5-lm) ground-based telescopes (Steffen et al. 2007). It 
is difficult for these telescopes to detect such low-mass planets by other methods, except via 
the microlensing technique (e.g., Yee et al. 2009). Such small telescopes are relatively plentiful 
and widely distributed around the world and thus their availability for observing phenomena 
such as multiple transit timings is more flexible than that for larger aperture telescopes. 

Searching for planets in MMRs with hot J up iters is also important. According to the 
core accretion models, hot Jupiters are predicted to form at ~ several AU where solid material 
is abundant enough to become a massive core and to accrete gas, and then migrated inward 
to their current positions in some way; mainly either via gravitational disk-planet interaction 
models (hereafter disk-planet interaction models, e.g., Lin & Papaloizou 1986), or via planet- 
planet scattering and/or Kozai process followed by tidal evolution (hereafter planet-planet 
scattering models, e.g., Nagasawa et al. 2008). The disk-planet interaction models predict that 
two planets including a gas-giant planet in a system will be easily captured in low-order MMRs 
during their migration processes (e.g., Thommes 2005). On the other hand, the planet-planet 
scattering models are thought to be difficult to form planets captured in MMRs. The current 
distribution of the sky-projected angle between the stellar spin and planetary orbital axes 
measured via the Rossiter-McLaughline effect shows that a significant amount of hot Jupiters 
are misaligned (e.g., Narita et al. 2009; Winn et al. 2009b; Triaud et al. 2010, hereafter T10). 
This fact implies that the planet-planet scattering models might play an important role in 
forming hot Jupiters. Furthermore, Winn et al. (2010) pointed out that hot Jupiters around 
hot stars tend to have high obliquities rather than those around cool stars, and they proposed 
a hypothesis that most or all hot Jupiters initially have high obliquities in some way and only 
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cool stars have damped their obliquities. If this explanation is true, planets in MMRs with hot 
Jupiters would be rare. The TTV method can thus be a powerful tool to check the existence 
or nonexistence of such planets around hot Jupiters and to test these planetary migration 
scenarios. 

Although the TTV method has been employed in many searches for additional planets, 
to date most of them have not shown any planetary signals while they have been used to put 
upper limits on masses of hypothetical additional planets in the systems (e.g., Adams et al. 
2010, and a summary therein). Recently, three systems have been reported to show plausible 
TTV signals; WASP-3 (Maciejewski et al. 2010a), WASP-10 (Maciejewski et al. 2010b), and 
Kepler-9 (Holman et al. 2010). 

The transiting planet WASP-5b was discovered by Anderson et al. (2008) (hereafter A08) 
as a 1.58 Mj up hot Jupiter orbiting a V = 12.3 (G4V) star with a period of 1.63 days. Gillon et al. 
(2009) (hereafter G09) have conducted high-precision photometric and spectroscopic follow-up 
observations for this system using the 8.1 m Very Large Telescope (VLT). Although the resulting 
photometric data suffered from uncorrectable systematic errors and were not used for their 
analysis, they reanalyzed photometric and RV data presented in A08 and found marginal non- 
zero eccentricity (~ 2a) for the planetary orbit. They also found that the reduced-^ 2 of a linear 
fit for four transit timing data was 5.7, which might be a sign of additional bodies. Southworth 
et al. (2009) (hereafter S09) presented two high quality photometric transit light curves obtained 
with the 1.54 m Danish telescope and with residual standard deviations relative to theoretical 
fits of 0.50 mmag and 0.59 mmag. They argued that the inconsistency with a linear ephemeris 
could be associated with a relatively poor transit light curve from the Faulkes Telescope South 
(FTS). T10 presented 33 new RV data obtained with the HARPS spectrograph, including data 
during a transit, and 5 CORALIE data in addition to the 11 given by A08. They derived the 
sky-projected spin-orbit angle of A = 12.1°tij) o°, which is consistent with a spin-orbit alignment 
that is naturally expected from the disk-planet interaction models for the migration mechanism 
of WASP-5b. However, the planet-planet scattering models followed by the stellar obliquity 
dumping could still be an alternative scenario, because the host star is a "cool star" in the 
context of Winn et al. (2010), which might therefore easily dump the stellar obliquity. We 
cannot exclude this scenario at this point in time, and therefore searching for additional plants 
in low-order MMRs via the TTV method is useful to test these migration scenarios. 

In this paper, we present observations of 7 new transits of WASP-5b obtained with a 
61 cm telescope. We also gathered available photometric and RV data and performed joint- 
fit to refine transit parameters and to search for TTV signals. In section 2 we describe our 
observations and data reduction. We show our light curve modeling in section 3 and discuss 
results in section 4. We summarize our findings in section 5. 
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2. Observations and Data Reduction 



2.1. Transit Observations 

We observed 7 transits of WASP-5b (corresponding transit epochs, E, are: 160, 244, 
432, 451, 459, 607, and 615 based on the ephemeris given by A08) using the 61 cm Boiler 
and Chivens (B&C) telescope located at Mt. John University Observatory, operated by the 
University of Canterbury at Lake Tekapo in New Zealand. The B&C telescope is normally 
used for the microlensing followup program of Microlensing Observations in Astrophysics col- 
laboration (MOA; Bond et al. 2001; Sumi et al. 2003, 2010), and the transit observations were 
obtained during low-priority observing time. A lk x lk CCD camera of Apogee ALTA U47 is 
mounted on the telescope which has a focal length of 8.3 m. The CCD pixel size of 13 /im x 
13 fim therefor leads to a pixel scale of 0''33 per pixel and a field of view of 5 '5 x 5f5. The 
CCD frame readout time is 5 s and has a readout noise, ri read , of 11.5 electrons per pixel. The 
dark current noise, n dar k, is 0.03 ADU per pixel per second under typical conditions. A Bessell 
I filter was used for all observations and the exposure times were 30 s for £"=245 and £"=607, 
and 60 s for the others. All images were taken with the telescope properly focused, except for 
the transit of £=459 when the telescope was slightly defocused so that the FWHM of the PSF 
was 9-12 pixels (3."0-4"0). The typical seeing at the observatory was ~ 2."0. An observing log 
is shown in table 1. 

2.2. Data Reduction 

All images were initially bias subtracted and flat-field corrected in the standard manner. 
We then performed synthetic aperture photometry for the target star WASP-5 and some (2 - 
4) comparison stars in the same field of view, using the following procedure. Note that each 
star to be measured is enough isolated from the closest star not to be contaminated during the 
procedure. 

First, for each image j and using a step size of 0.1 pixels, we searched for an initial 
photometric aperture radius r opt j that maximized the signal to noise ratio of the target flux. 
The signal, -F s tar, is flux from the target which is equal to the total flux in the aperture minus the 
corresponding sky flux, F sky = m pix / s ky. Here, m pix is the total number of pixels in the aperture 
and / s ky is the median flux per pixel in an annulus between radii r opti j + 50 and r opti j + 60 where 
the flux contribution from the target star is negligible. The total noise contributing to the 
signal is modeled as: 

^total = \J A^tar + ^sky + ^read + ^dark + ^L) (!) 

where iV star = \J F sta _ r / g is the photon noise of the target star, iV s ky = sj ' F sky / 9 is the sky back- 
ground noise, iV read = y/m pix n reaid /g is the read out noise, iV dark = y/m pix n d£irk At is the dark 
current noise, and iV SC i n is the atmospheric scintillation, respectively. Here, g = 1.19 electron 
per ADU is the CCD gain and At is the exposure time in seconds. n read and n dar k are defined 
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in the previous section. The scintillation noise can be expressed as 

z 7/i ( h \ 

N "*» = tr °£,2/3 (At) l/2 eX P (~ 800to J Fstar ' (2) 

where z is the airmass, D is the telescope diameter in cm, and h = 1029 m is the observatory 
altitude (Dravins et al. 1998; Young 1967). <jq is a coefficient which is often taken to be 0.064 
and we also adopted this value. At typical conditions with exposure time of 60 s, the main 
contributor to the total noise was the photon noise (~0.0012% of the target flux) and the 
contributions from the other noises were one or two orders of magnitudes less than that. 

Next, after determination of the r opt j value for the target star, this radius was applied 
to all other comparison stars on the same image to measure their fluxes. The target flux was 
then normalized to a reference flux that was created as the weighted average of the fluxes of the 
comparison stars. The photometric error of the normalized flux was initially estimated using 
the equation (2) and the error propagation equation. 

Finally, for each transit E, we searched for a scaling factor 7^ (in the range 0.7 to 1.50) 
to the initial radius r op t,j in order to minimize the rms scatter of the derived out of transit 
(OOT) light curve (either before or after transit). 

In order to assess the validity of our method, we also performed aperture photometry 
with the fixed radius rg Xi £ for each transit which was selected so as to produce the minimum 
rms scatter for the OOT light curve. The result of this exercise was that the OOT rms values 
derived from the ^E r opt,j apertures were slightly better than those from the apertures 
for all transit light curves, except for the transit £'=615 for which there are a relatively small 
number of OOT data points, possibly leading to a large statistical fluctuation. Table 2 shows the 
comparison of OOT rms values derived from the two methods. In order to maintain consistency, 
all light curves were generated by using the scaled apertures method and these were used for 
further analyses. 

All time stamps of observations, which were recorded in the FITS headers as the obser- 
vation start and end points in units of Julian Day (JD) based on Coordinated Universal Time 
(UTC) and obtained from the network time protocol, were converted to the mid-observation 
points in units of Barycentric Julian Day (BJD) based on Barycentric Dynamical Time (TDB) 
using the code UTC2BJD (Eastman et al. 2010). The time given the network time protocol 
has occasionally been confirmed to coincide with the GPS time within 1 s by observers. The 
time-conversion error does not exceed 20 ms. 

3. Light Curve Modeling 

3.1. Systematic correction 

Since apparent systematic trends were seen in the derived light curves, we corrected these 
systematic effects as follows. The systematic trends could arise from the changing airmass, 
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variations in the atmospheric extinction coefficient, slow brightness variability of the target or 
comparison stars, and so on. Using the correction factor Am corr in the magnitude scale, the 
corrected flux £ corr can be expressed as 

Fcorr = Fobs X 1(T - 4A — (3) 

where £ OD s is the observed flux of the target normalized by the reference flux. If we assume that 
amplitudes of the intrinsic stellar variability and the change in the second-order atmospheric 
extinction coefficient are both proportional to time, Am corr can be expressed as 

Afflcon = k + k z z + k t t + k tz tz (4) 

where z is the airmass, t is the time, and (ko, k z , k t , k tz ) are the relevant coefficients. This 
equation is similar to equation (1) of Winn et al. (2009a), although here we use four coefficients 
instead of three (see Appendix). Provided there were enough data points, only the OOT light 
curve segments were fitted to estimate the systematic-correcting function. However, three of 
the light curves (£=244, 432, and 607) did not have an adequate amount of data points in their 
OOT sections, and also apparent systematics were seen in the first half part of the light curve 
of £'=451 transit which could not be corrected if only the OOT data were used. Therefore, we 
decided to fit each light curve including the transit section while fixing the parameters defining 
the transit shape to the values derived by S09, who had obtained two transit light curves with 
higher precision than ours. 

To create the parameterized transit light curve, an analytic model of Mandel & Agol 
(2002) was used. This model requires 6 parameters: the planet/star radius ratio R p /R*, the 
planetary orbit semi-major axis in units of the star radius a/R*, the orbital inclination to the 
line of sight i, the orbital period P, and two stellar limb-darkening coefficients U\ and u-i- We 
used a quadratic limb darkening law, 

J(/i) = l-n 1 (l- /U )-n 2 (l- /U ) 2 , (5) 

where I is the intensity and fi is the cosine of the angle between the line of sight and the 
line from the center of the star to a position of the stellar surface. For correcting systematics, 
these parameters were fixed to the values of R p /R* = 0.1110, a/R* = 5.4142, % = 85.8°, P = 
1.6284246 d (adopted from S09), u 1= 0.257, and u 2 =0.341 (adopted from the tables of Claret 
2000). Using this model and the equation (3) and (4), each light curve was fitted using the 
X 2 minimization method with five free parameters: ko,k z ,k t ,k tz , and the mid-transit time T C £. 
The corrected light curves are shown in figure 1. 

To check for the validity of this correction, we alternatively corrected the £=160 com- 
plete transit light curve by fitting only its OOT parts, then compared the two corrected light 
curves. Figure 2 shows the two light curves and their differences. The difference was no larger 
than 0.00048, which is well within the typical error bar that was estimated from equation (2) 
(~0.002). Since we obtained similar results for two other complete transits (£=459 and 615), 

7 



we decided to use the light curves corrected by the former method for further analyses. 

Next, to estimate a statistical weight for each light curve, we fitted each one separately 
and rescaled the initial error bars so that each reduced-^ 2 became unity. We also estimated 
time-correlated systematic noise (red noise) by using "time-averaged" method according to 
Winn et al. (2008). In this method, an error scaling factor /3 ( > 1) is calculated for each light 
curve by analyzing the residuals from the best-fit model. If the red noise is not significant, the 
standard deviation after binning the residuals into M bins of N points, a^, would be expected 
to be 



where (j\ is the standard deviation of the un-binned residuals. The actual measured value of 
o"at is usually larger by a factor (3. Since the value determined for (3 depends on the choice of 
N and M, we adopted the maximum /3 value in the range from iV = 5 to 30 and rescaled each 
error bar by this value. The derived values of f3 and entire rms scatters of residuals are shown 
in table 3. 

3.2. Joint Fit 

To assess the quality of our light curves compared to the previous work, we fitted the 
seven transit light curves from the B&C telescope simultaneously and compared the resultant 
parameters to the values determined by S09 and T10, in which independent data sets were used. 
At this time we used the same parameterization as previous subsection and used formulae of 
Ohta et al. (2009) for creating a transit light curve model. We used common parameters of 
i, Rp/R*, a/R*,Ux, and w 2 to all transit light curves, while used respective transit-mid times, 
T c ,e- We allowed the parameters of i, R p /R*, a/R*, u\, and each T c< e to vary, while u<i was 
fixed at the theoretical value of 0.321. We fixed U2 in the fitting process since the u\ and u<i 
quantities are strongly correlated and could not both be well determined from the light curve 
fitting. The eccentricity was also fixed to zero. 

We then derived the best fit parameters by minimizing the \ 2 statistics using the AMOEBA 
algorithm (Press et al. 1992) and estimated uncertainties using the A% 2 = 1.0 criterion, following 
Narita et al. (2007). The resultant parameters and their l-a uncertainties are shown in table 
5. The parameters from the B&C are consistent with the ones derived by S09 and T10 within 
their error bars, except for the period P, which may be a sign of the TTVs (discussed in section 



Figure 3 compares phase folded, 120-s binned light curves from the B&C and the Danish 
(online data of S09) telescopes. Before folding the Danish data, the quoted errors were normal- 
ized and rescaled using the same method for the B&C data descried in the previous section, 
i.e., normalized so that the reduced-^ 2 for the best- fit to each transit was unity and rescaled 
by a factor /? accounting for the red noise; we then took the weighted average and its error for 




(6) 



4). 



8 



each binned data set. 

The rms residuals from the best-fit models for the Danish data and the B&C data are 
0.00046 and 0.00081 respectively. However, uncertainties in the derived transit-model param- 
eters from the B&C data are comparable to the ones from the Danish data presented in S09 
(see table 5). This result initially appears surprising given the difference in the fit residuals; 
however, we found that the values of the error-rescaling factor for the Danish data are larger 
than the most of that for the B&C data (see table 3). As a consequence, the average values of 
the binned errors of the Danish and B&C data sets (0.00061 and 0.00084 respectively) are more 
comparable. Moreover, the i?-band filter used for the Danish observations would have resulted 
in a larger limb-darkening effect than that seen using the J-band filter employed in the B&C 
observations. This leads to a relatively poorer determination of the model parameters for the 
Danish data (e.g., Pal 2008). In addition, S09 included limb-darkening model dependencies in 
their error estimations which enlarged the errors of the transit-model parameters. However, all 
general limb- darkening models produce a symmetrical transit shape, which has little effect on 
errors of transit timings. For this reason, we do not include model dependencies derived from 
limb-darkening in this paper. 

To refine transit model parameters and to revise the timings of the previous transits, 
we assembled all available photometric and RV data in addition to our own data, and fitted 
them jointly. The photometric dataset consists of 11 transits including one transit from the 
1.2 m Euler telescope (133 data points) and one transit from the 2.0 m FTS telescope (335 
data points) analyzed in A08 (private communication), two transits from the 1.54 m Danish 
telescope (174 data points) presented in S09, and seven transits from the B&C telescope (1870 
data points). The RV dataset consists of 16 and 33 data points from CORALIE and HARPS 
instruments, respectively, presented in T10. In these available datasets, all time stamps are 
provided either in the form Heliocentric JD (HJD) (Euler, FTS, and Danish) or BJD (CORALIE 
and HARPS). However, time-standards (e.g. UTC or TDB) on which these time systems are 
based are unspecified in the publications. As Eastman et al. (2010) alerted recently, specifying 
the time-standard is important in order to achieve 1 minute accuracy, and TDB-based BJD 
is the recommended time system. We confirmed (in private communications) that all time 
systems in the available data were based on UTC , and we therefore converted all of them to 
the TDB-based BJD. 

In order to treat the acquired photometric data equally with the B&C data, their quoted 
errors were normalized and rescaled in the same manner as for the B&C as explained in the 
previous subsection. The derived value of (3 and the rms scatter of the residuals for each transit 
are shown in table 3. For the RV data, systematic errors of ~7 m s _1 may be presented (A08 
and T10) in addition to the quoted internal errors, these often come from stellar activity (stellar 
jitter). To add statistical weight to the RV data, we rescaled the quoted RV errors for both 
CORALIE and HARPS by adding 7 m s -1 quadratically. This value is consistent with the 
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upper value from an empirical model of Wright (2005) for a G-type star. 

We then fitted these data with three additional parameters: the RV semi-amplitude K, 
the sky-projected stellar rotational velocity ^sin/ s , and the sky-projected spin-orbit alignment 
angle A. The last two parameters were needed to model the Rossiter-McLaughline(RM) effect 
during transits of WASP-5b. The RM formula we used was based on Hirano et al. (2010), which 
is more appropriate than the old one used in the previous work (T10). We allowed the param- 
eters of i, R p /R*, a/R*, T Ct E, K, ^sin/ s , and A to be free. We also fitted u\ for each telescope 
while each U2 was fixed at the theoretical value (Claret 2000, 2004), because, as mentioned 
previously, it was difficult to determine unique values for both u\ and u<i simultaneously due 
to their strong correlation. As for the eccentricity, we first allowed it to be free and obtained 
the value of 0.002 

iaoo2- This value is consistent with zero within 1-a and also consistent with 
the result of T10 (e < 0.0371, 2 a), but marginally inconsistent with G09 (e = 0.049j£g??). 
However, the latter used a smaller RV data set in comparison with us, and hence we adopted 
zero for the eccentricity and re-fitted the data. The derived value of U\ and the fixed value of 
U2 for each telescope are listed in table 4, and the other best-fit parameters are shown in table 
5. 

From the combined fit, we refined some transit model parameters, i, R p /R*, and a/R* 
against the published values while they are consistent within their error bars. Uncertainties in 
K and Vsml s became larger than the ones determined by T10 because we incorporated the 
stellar jitter in the RV errors. We adopted the value of the stellar mass as the one determined 
spectroscopically by T10, and used this to convert fitted parameters to physical parameters. 
The derived values are shown in table 7. The x 2 values for the COLALIE and HARPS datasets 
are 16.0 and 32.7 respectively, which are very close to the number of data points (16 and 33 
respectively). 

We confirmed that the projected spin-orbital angle A is consistent with zero (A = 
7.2°l;g;5o), which was first reported by T10, by adding independent photometric data (the 
B&C and Danish data) and using an improved RM formula. 

4. Results and Discussion 

4-1. Significance of TTV signal 

We show the transit timings with 1-a uncertainties in table 6. Using these data, a new 
ephemeris was calculated via a linear fit using a function 

T C , £ = T C , + £P, (7) 

where T c , is a reference time of E=0. The results are T c>0 [BJD TD b] = 2454375.62510 ± 0.00019 
and P = 1.62843142 ± 0.00000064 d. This fit yields x 2 = 32.2 for 9 degrees of freedom which 
indicates that a linear fit does no mach at a confidence level of 99.982 %, or 3.7 a. Figure 
4 shows the timing residuals and their uncertainties. This result implies that we see excess 
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variations in the transit times due either to unknown systematic effects or possibly to real 
TTVs. A liner fit for only the B&C data yields x 2 = 17. 1 for 5 degrees of freedom which 
corresponds to a confidence level of 99.57 % or 2.9 a, while the fit for the other 4 data yields x 2 
= 9.8 which corresponds to a confidence level of 99.27 % or 2.7 a. Even if two data having the 
largest uncertainties in all 11 timing data (corresponding epochs of E=7 and 451), which may 
have larger systematics, are discarded, the x 2 value for a linear fit remains 27.1 for 7 degrees 
of freedom, which corresponds to a confidence level of 99.968 % or 3.6 a. This result increases 
the statistical significance of a TTV excess implied by G09, who used four timing data (two of 
them were derived from same photometric data as we use) and derived X/dof f° r a linear fit of 
5.7 which corresponds to a significance of 99.67 % or 2.9 a. 

The standard deviation of the observed timing residuals is 68 s and the mean uncertainty 
of the timings is 41 s, therefore if this is a real TTV signal, the actual standard deviation of 
the TTV is expected to be ~50 s. Such a large timing deviation could not be explained by 
effects other than additional perturbing bodies. The size of these alternative possibilities are: 
the Applegate effect (~ 1.5 s per 11 years, Watson & Marsh 2010), the light travel time effect 
due to an outer massive body 1 s per 3 years, e.g., Agol et al. 2005), orbital decay caused 
by tidal dissipation (~ 5 ms year -1 , Hellier et al. 2009), and orbital precessions due to tidal 
deformations (< 130 ms year -1 ), due to the general relativity effect (< 2 ms year -1 ), and due 
to stellar quadrupole moment (< 0.005 ms year -1 , Heyl Sz Gladman 2007; Jordan & Bakos 
2008) (see table 8). In addition, the presence of an exomoon also could not be responsible for 
the signal, because Weidner & Home (2010) placed an upper mass and distance (9.1 xlO -4 M® 
and 4.39 Rj up respectively) of a possible moon for the WASP-5b system by considering the 
three-body stability issue; such moon would result in only ~0.02 s TTVs. Thus, the remaining 
possibilities for explaining the signal are additional planets or Trojan companions (bodies at 
1:1 MMR). 

4-2. Upper Mass Limits for Secondary Planet 

A marginal evidence for TTV signal is seen, however, it is hard to determine the char- 
acteristics of the additional body (planet) at this time because of a large parameter space for 
the planet parameters: mass, orbital period, eccentricity, position of periapsis, reference phase, 
and mutual inclination. Alternatively, using the observed TTV data, we place upper limits on 
the mass of the hypothetical additional planet as a function of period ratio of the two planets 
(WASP-5b and the secondary body) by numerical thee-body simulations, similar to several 
previous analyses (e.g., Steffen & Agol 2005; Miller-Ricci et al. 2008a,b; Bean 2009; Gibson 
et al. 2009; Adams et al. 2010; Hrudkova et al. 2010; Gibson et al. 2010). 

In order to simplify things, we assume that the secondary planet orbits in the same 
plane as WASP-5b (the primary planet) and that the two planets have initially circler orbits. 
Because the amplitude of timing variations tends to increase as the orbital eccentricity of the 
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second planet increases, basically we only need to consider an initially circular orbit for the 
secondary planet (e.g., Gibson et al. 2010). For a given initial orbital period, mass (M2), and 
initial phase (82) of the secondary planet, we perform numerical integrations of the equations 
of motion for the three-body system using the 4th-order Runge-Kutta method and a time step 
of 1 s. We continue the integrations for an equivalent elapsed time of 3,000 days, which is three 
times longer than the observation time interval. When the primary planet in the numerical 
procedure passes through the star-observer field of view the time step is reduced to 0.001 s and 
the calculations are used to create simulated transit timing data. These data are fitted to a 
linear function in order to obtain an "observed" orbital period and synthetic TTV data, which 
are then compared to the observed TTV data and a \ 2 value obtained. Transit timings of the 
secondary planet are also extracted in order to calculate its observed orbital period. 

For a given secondary-planet model (initial period ratio, mass, and phase), a x 2 value 
is chosen so as to be a minimum by shifting the reference epoch of the artificial TTV data 
sequentially. We also assume that the TTV amplitude is proportional to M 2 (Agol et al. 2005; 
Holman & Murray 2005) and integrate only the case of Mi = 3 x 10~ 6 M® for each initial period 
ratio (ranging from 0.2 to 5.0, increased by a factor 1.005) and 9 2 (stepped by 30 degree), then 
rescale M 2 (or the amplitudes of the artificial TTV data) so that the A% 2 , between the x 2 value 
derived above and a \ 2 value for linear fit, becomes 9.0 (3 a confidence limits). 

To check the long term stability, we also conduct long term runs for up to 10 5 years with 
a time step of 100 s and eliminate models in the cases the system becomes unstable (one body 
is ejected or two bodies collide). We then derive 12 upper-limit values of M2 from 12 models of 
62 as a function of observed period ratio, and derive a maximum upper-limit value by taking a 
maximum value for a given period ratio from the linearly-interpolated upper-limit lines. 

The derived 3 a upper limits are shown in figure 5 (black solid line). The overplot- 
ted thick-dashed line represents the boundary of Hill-stable regions calculated from Barnes 
& Greenberg (2006), inside which long-term stabilities are not guaranteed while some stable 
regions can still exist (Barnes & Greenberg 2006). We also overplot an upper-limit line esti- 
mated from the RV data (dotted line), which corresponds to a line causing the RV amplitude of 
21 m s _1 ; we assume that the added systematic error of 7 m s _1 is entirely due to the possible 
additional body and multiply this value by three as an upper limit. As a result, outside the 
boundaries, we place upper limits on the secondary mass as 22-70 M® in the period ratio from 
0.2 to 5.0 from the RV data, and more stringent limits down to 2 M® near 1:2 and 2:1 MMRs 
from the TTV data (at the 3 a level respectively). 

We also put the upper limit on the mass of a Trojan companion (an object at one of 
the two triangular Lagrange points, L4 and L5, of the planet-star orbit) as 100 M® (3 a) from 
the TTV data. The upper limit of ~40 M® at 1:1 MMR from the RV data is not correct, 
because if a Trojan companion exists at L4 or L5 points, the observed RV data would show one 
planet orbiting at the barycenter of the primary planet and the Trojan companion. Instead, 
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a method to find an imbalance of mass at the L4/L5 points using both RV and photometric 
data has been proposed by Ford & Gaudi (2006). The basic idea of this method is to observe 
a difference between the time of vanishing stellar RV variation (T^kv) an d the time of the 
mid-transit (T c ). Madhusudhan & Winn (2009) estimated the upper mass of the imbalance 
mass of Trojan companions in the WASP-5 system as <54.7 M ffi (2 a) at the L5 point (behind 
the planet), based on the data of A08. Here we improve this estimate by employing more RV 
and photometric data. We measure the time difference At = T C) o - T^ry with (without) RV 
data during transit (14 data points) adopting the eccentricity of WASP-5b is zero, resulting At 
= 4.8 +/- 4.4 min (12.8 +/- 4.8 min). Accordingly, we set upper limits on the excess mass of 
the Trojan companions, Mr, which is defined as the difference in the mass at L4 (Mr,L4) and 
the mass at L5 (Mt,l,5) (namely, Mr = Mt,l,4 — Mt^), through the relation; 



Madhusudhan & Winn (2009) and originally from equation (1) in Ford & Gaudi (2006) ]. We 
find Mp = 7.4 ± 6.8 (19.8 ± 7.7) M ffi and put an upper limit on the excess mass near the L4 
point as ~28 M ffi (~43 M e ) at the 3 a level. This result lowers the limit derived from only the 
TTV data (< 100 M e ), and improve the one derived by Madhusudhan & Winn (2009) by at 
least a factor 2. 

4-3. Example Models of Secondary Planet 

Here, in order to illustrate example models of secondary planet which can account for the 
observed timing variations, we fit the simulated TTV data to the observed one and search for 
well-fitting models. Using the artificial TTV data generated in the previous section, in which 
the initial eccentricity of the secondary planet, e^-, was set to zero, we search for a best-fit 
solution for each orbital period ratio by scaling TTV amplitude and shifting a reference epoch 
of the artificial TTV data. At this time, we use only the case that the initial phase of the 
secondary planet is zero for simplicity. Since, generally, a libration period of TTVs increases 
as the orbital period ratio approaches to a MMR, solutions around MMRs should be searched 
finely. For this reason, we generate a set of additional artificial TTV data around the low-order 
MMRs of 1:4, 1:3, 1:2, 3:5, 5:3, 2:1, 3:1, and 4:1, by 10 times denser than the other regions, and 
search for best-fit solutions also for the additional period ratios. We then convert the best-fit 
TTV amplitudes to the secondary masses, assuming a TTV amplitude is proportional to a 
secondary-planet mass. In addition to the case of e2=0, we also generate a set of artificial TTV 
data for the case of e2=0.1 in the same conditions as the previous one and in an additional 
condition; initial phase of the periapsis of the secondary planet which is set to zero. In this 
case, the regions around the MMRs of 1:4, 1:3, 2:5, 1:2, 3:5, 2:1, 5:2, 3:1, and 4:1 are searched 
densely. We then search for the best-fit solutions against the respective orbital period ratios. 




(8) 



where P and M p are the orbital period and mass of WASP-5b respectively [ equation (2) in 
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Figure 6 shows the derived best-fit mass of the secondary planet as a function of the orbital 
period ratio for the cases of e2=0 (left top panel) and e2=0.1 (right top panel), and their x 2 
values (bottom panels). 

As a result, for the case of e2=0, the regions where the best-fit masses are lower than the 
upper limits placed from the radial velocity data (RV limits) are limited only around the low- 
order MMRs of 1:3, 1:2, 3:5, 5:3, 2:1, and 3:1. Such regions then extend to around other MMRs 
as C2 increases to 0.1. In addition, the \ 2 maps show there are a number of local minima over 
the period-ratio range, reflecting the fact that the libration period of artificial TTVs gradually 
changes as the orbital period ratio changes, and some of the local minima around MMRs have 
low-x 2 values, less than the degrees of freedom of 9. Figure 7 shows zoom around the 2:1 MMR 
of figure 6 for example. These facts suggest that it is possible to explain the observed TTV 
data by a certain model of a perturbing planet having lower mass than the RV limit, and such 
a model would likely be near low-order MMRs, if its eccentricity is enough small. Non- or 
a small eccentricity of the secondary planet is naturally expected if the planet has migrated 
by the disk-planet interactions, and if the above scenario is true, this fact could be important 
observational evidence that the hot Jupiter has been formed according to one of the disk-planet 
interaction models (e.g., Thommes 2005). This would also be consistent with the fact that the 
sky-projected spin-orbital angle is consistent with zero. 

We show four example TTV models (a, b, c, and d), having lower masses than the 
RV limits and locating near MMRs, with the observed TTVs in figure 8. The locations in 
parameter space of period ratio and mass are marked as star symbols in figure 6. (a) is the 
least x 2 model (x 2 =4.35) in all models considered here, which is located near the 1:2 MMR 
and 0.003 Mj up in the case of e2=0. (b) is a low-% 2 model (x 2 =5.01) located near the 2:1 
MMR and 0.093 Mj up in the case of e2=0. Since (a) and (b) show similar libration periods and 
amplitudes, distinguishing them by only ground-based TTV observations would be difficult 
and radial velocity followup observations would be necessary, if one of such models is true, (c) 
is another low-% 2 model (x 2 =6.65) located near the 2:1 MMR and 0.011 Mj up in the case of 
e2=0.1. (d) is a model located near the 4:1 MMR and 0.053 Mj up in the case of e2=0.1 with 
X 2 =20.57. The \ 2 value of the model (d) is somewhat high due to the relatively longer libration 
period, however, if two timing data having the largest error bars (E=7 and 451) are discarded, 
the x 2 value decrease to 11.6 for 7 degrees of freedom, and therefore such longer libration-period 
models are also thinkable. 

On the other hand, non-MMR with rather higher eccentricity models might also be 
possible for which we do not search here. If this case is true, it would be the first example 
where a hot- Jupiter system has a non-MMR, low- mass planet with a short orbit. Recently it 
has been revealed that most multi-planetary systems consisting of several hot super-Earths are 
not in MMRs (Mayor et al. 2009a,b) and this can not be explained naturally using standard 
migration models (Terquem & Papaloizou 2007; Ogihara et al. 2010). Thus, it is important to 
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check whether the situation in hot- Jupiter systems is same or not. 

In order to distinguish these models and find an exclusive solution, additional high- 
precision and extensive timing data are necessary, and high-precision radial velocity followup 
observations are also helpful to confirm the additional planetary signal and constrain the phys- 
ical parameters. A successful program would provide valuable information for improving our 
understandings of planetary formation scenarios. 

5. Summary 

We have observed 7 new transits of the hot Jupiter WASP-5b using 61 cm telescope 
in order to search for an additional body via the TTVs method. By combining all available 
photometric and RV data, we slightly refined the transit-model parameters and confirmed all 
parameters were consistent with published values, except for the orbital period which may be 
a sign of TTVs. We also confirmed that the sky-projected spin-orbit angle is consistent with 
zero, which was first reported by Triaud et al. (2010), by including independent photometric 
data and using an improved RM formula. This fact is consistent with a disk-planet interaction 
model being the favored migration scenario for WASP-5b, although planet-planet scattering 
models cannot be excluded at this point in time. 

The x 2 value for a linear fit to the 11 transit timings is 32.2 for 9 degrees of freedom. 
This result indicates that the transit timings do not match a linear ephemeris at a confidence 
level of 99.982 %, or 3.7 a, due either to unknown systematic effects or possibly to real TTVs. 
If this signal is real, the standard deviation of the TTVs is as large as 50 s, and the most likely 
cause of this is the presence of an additional body in the system. From the RV data, we put an 
upper limit on the RV amplitude caused by the possible secondary body (planet) as 21 m s" 1 , 
which corresponds to its mass of 22-70 M ffi over the orbital period ratio of the two planets 
from 0.2 to 5.0. From the TTV data, using the numerical simulations, we place more stringent 
limits down to 2 M® near the 1:2 and 2:1 MMRs with WASP-5b at the 3 a level, assuming the 
two planets are co-planer. We also put the upper limit on excess Trojan mass as 43 M e (3 a) 
using both RV and photometric data. 

We also find that if the possible secondary planet has non- or a small eccentricity, it would 
likely be near low-order MMRs. The fact would also support the disk-planet interaction models 
as the migration mechanisms for the hot Jupiter. Alternatively, if the secondary body is an 
eccentric, non-MMR planet, it would be a challenge for the current planetary migration models. 
Further followup observations for the WASP-5 system by both photometry and spectroscopy 
will reveal the true nature of the TTV signal and shed light on the migration mechanisms of 
planetary systems. 
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Appendix. 

When we consider primary and secondary atmospheric extinctions, time-averaged stan- 
dard magnitudes of a target star and a reference star in an arbitrary passband, M ob j and M re f, 
can be expressed as 

M obj = m obj - k'z - k"C oh] z + TC bj + Z P + / obj (t) , ( Al ) 

M ref = m ref - k'z - k"C rcl z + TCW + Z v + / ref (t) , ( A2) 

where m ob j and m rcf are observed magnitudes of the target star and the reference star, C ob j 
and C rc f are their color indices, z is the airmass, T is the transformation factor, Z p is nightly 
zero point, k! and k" are coefficients of primary and secondary extinctions, and f h s {t) and 
f re i(t) are terms of stellar intrinsic variability as a function of time (see e.g., Chapter 6 of 
Warner 2006). When we approximate that each stellar intrinsic brightness varies with time in 
linear function, then the differential magnitude between the target star and the reference star, 
AM = M obs — M re f , is written as 

AM = Am - k"ACz + TAC + at + b, (A3) 

where Am = m ob j — m rc f , AC = C ob j — C rc f , and a and b are coefficients. Here, we assume that 
the k" varies with time t in linear function, as k" = k" + k'\t, then 

AM = Am - (k" + k'\t)ACz + TAC + at + b. (A4) 

Consequently, if we redefine k z = —k" AC, k tz = —k'^AC, k t = a, and k = TAC + b, then the 
correction magnitude Am corr = AM — Am can be written as 

Am corr = k z z + k tz tz + k t t + k . (A5) 
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Table 1. The observing Log 



Date 


epoch (E) 


exp. [sec] 


niter 


fp oi data 


iocus/detocus 


airmass 






2008 June 18 


160 


60 


I 


248 


tocus 


1.55 -> 1.00 




zUUo fNovember z 


Z44 


oU 


1 




focus 


1.00 


1.35 




2009 September 4 


432 


60 


I 


168 


focus 


1.03 


1.00 - 


> 1.07 


2009 October 5 


451 


60 


I 


313 


focus 


1.07 -> 


1.00 - 


> 1.39 


2009 October 18 


459 


60 


I 


248 


dcfocus 


1.01 


1.80 




2010 June 16 


607 


30 


I 


307 


focus 


2.60 


1.23 




2010 June 29 


615 


60 


I 


202 


focus 


1.78 


1.04 





Table 2. Comparisons of rms scatters between OOT light curves derived using optimized radii (7_Er opt j) and that derived 
using fixed radii (rfi x E)- Noot is the number of OOT (either before or after the transit) data points and "b" or "a" in 
parentheses stands for before or after the transit, respectively. 



E 


Noot 


IE 


lET ov t,j 




rms scatter for OOT (%) 








(pixels) 


(pixels) 


lEr op t,j 


Tfix.E 


160 


75 (a) 


0.96 


9.6 - 13.1 


13.1 


0.177 


0.178 


244 


119 (a) 


1.27 


13.1 - 18.5 


15.8 


0.488 


0.492 


432 


50 (b) 


1.34 


11.5 - 14.5 


14.1 


0.205 


0.205 


451 


77 (a) 


1.47 


14.4 - 16.0 


14.0 


0.282 


0.284 


459 


52 (b) 


1.24 


13.9 - 17.1 


13.9 


0.179 


0.181 


607 


49 (a) 


1.42 


12.5 - 14.3 


13.0 


0.287 


0.288 


615 


44 (a) 


1.30 


11.1 - 13.7 


12.7 


0.161 


0.158 
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Table 3. Estimated red-noise factor /}, the number of data points in one bin (N ma x), un-binned rms scatter from the best-fit 
model (<Ti), and the total number of data points, for each transit light curve are listed. 



E 


8 


1 v max 


Ch [%1 


=# of data 




160 


1.15 


15 


0.197 


248 


B&C 


244 


1.16 


14 


0.364 


386 


B&C 


432 


1.09 


6 


0.205 


168 


B&C 


451 


1.76 


18 


0.234 


313 


B&C 


459 


1.01 


8 


0.220 


247 


B&C 


607 


1.31 


17 


0.322 


306 


B&C 


615 


1.19 


13 


0.210 


202 


B&C 


5 


1.16 


10 


0.094 


138 


Euler 


7 


2.28 


27 


0.262 


335 


FTS 


204 


1.28 


15 


0.046 


73 


Danish 


218 


1.34 


15 


0.054 


101 


Danish 



Table 4. The quadratic limb-darkening coefficients (u\ and U2) for each telescope used for the final joint fit are listed. Each 
ui was allowed to be free while each 112 was fixed at the theoretical value extracted from tables of Claret (2000) or Claret 
(2004). 



Telescope 


filter 


u\ (fitted) 


U2 (fixed) 


Euler 


R 


n S3 + 05 

u -°° -0.05 


0.32 


FTS 


SDSS i" 


49 +0 12 

u -^ y -0.13 


0.32 


B&C 


/ 


n 07 +0.05 

u - z ' -0.05 


0.32 


Danish 


R 


n oa +0 04 
U ' J4 -0.04 


0.32 
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Table 5. The best-fit parameters and their 1-a uncertainties for WASP-5b derived from the joint fits for only the 7 light curves from the B&C telescope (the second 
column) and for all available RV and photometric data (the third column) are listed. Parameters derived by S09 (the fourth column) and T10 (the rightmost column) are 
shown for comparison. 



This work (B&C only) 



This work (all) 



Southworth et al. (2009) Triaud et al. (2010) 



P [days] 
To [BJDtdb -2450000] 

a/i?* 
K [m s _1 ] 
e 

VsmI [km s _1 ] 
A H 



1.6284301 ± 0.0000012 
4375.62589 ± 0.00052 



85.01 
0.1116 
5.26; 



+ 1.17 
-0.78 

f0.0010 
-0.0014 

-0.23 
-0.17 



(adopted) 



1.62843142 ± 0.00000064 
4375.62510 ± 0.00019 
85.58l ; 76 
0.1108 ± 0.0011 
5.37±0.15 
269.4 ±3.3 
(adopted) 
3.05 ±0.41 
7.2 ±9.5 



1.6284246 ± 0.0000013 
4375.62569 ± 0.00024 
85.8 ± 1.1 
0.1110 ± 0.0014 



5.41 



+0.17 
-0.18 



(adopted) 



1.6284229 



-0.0000044 
-0.0000039 



4373.996764 



+0.00015 
-0.00014 



86.2 



+0.8 
-1.7 



n 1 1 n^ +0 0019 

u.iiuo_ 0007 
5.49lg:g 
268.7lig 
< 0.0371 



3.24 



+0.34 
-0.35 



1 9 4+ 8 - 2 
iz -^-11.9 



Table 6. Transit Timings and their uncertainties. 





Transit Timing 
[UJUtdb - 2 450 000J 


1-c uncertainty 


telescope 


5 


4383.76751 


0.00028 


hauler 


i 


4387.02286 


r\ AAAO/ 1 

0.00086 


FTS 


160 


4636.17465 


0.00047 


B&C 


204 


4707.82531 


f\ AAAA1 

0.00021 


Danish 


218 


4730.62252 


0.00022 


Danish 


1 A A 

244 


4772.96212 


A A A A C 1 
0.00051 


B&C 


432 


5079.10849 


0.00044 


B&C 


451 


5110.04645 


0.00073 


B&C 


459 


5123.07627 


0.00041 


B&C 


607 


5364.08262 


0.00057 


B&C 


615 


5377.10969 


0.00048 


B&C 



Table 7. Derived physical values and their uncertainties for the WASP-5b system. The value of the stellar mass was adopted 
to the value derived by T10. The values presented in S09 and T10 are shown for comparison. 





This work 


Southworth ct al. (2009) 


Triaud et al. (2010) 


a [AU] 


0.02702 ± 0.00059 


0.02729 ± 0.00049 ± 0.00027 


A n97nn+0.00056 
u.uz 1 u»_ oo62 


R* [R©] 


1.082 ± 0.038 


1.084 ± 0.040 ± 0.011 


i n5fi+ 080 

l.UOO_o. 29 


M* [Mq] 


1.000 ± 0.065 


1.021 ± 0.055 ± 0.030 


l-000±g;gg 


P* [pe>] 


0.79 ± 0.10 


0.803 ± 0.080 ± 0.000 


o.84±S:S 


Rp [Rjup] 


1.167 ± 0.043 


1.171 ± 0.056 ± 0.012 


1 14+ - 10 

i - i4 -0.04 


M P [Mj„p] 


1.568 ± 0.071 


1.637 ± 0.075 ± 0.033 


1 ccc+0.067 

i.ooo_ 070 


Pp [Pjup] 


0.92 ± 0.11 


1.02 ± 0.14 ± 0.01 


1.05 ± 0.20 
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Table 8. Expected TTVs for WASP-5b induced by effects other than that of an additional perturbing planet or Trojan. 
T mod denotes the modulation time-scale of the host star. 



effect 


expected TTVs 


reference 


the Appregate effect 


— 1.5 (T morf /llyears) [s] 


Watson & Marsh (2010) 


light travel time effect 


< 1 [s/3ycars] a 


e.g., Agol ct al. (2005) 


exo-moon 


< 20 [ms] 


Kipping (2009); Weidner & Home (2010) 


tidal dissipation 


~ 5 [ms year -1 ] b 


Hellier et al. (2009) 


orbital precessions due to ... 






tidal deformations 


< 130 [ms year -1 ] c 


Heyl & Gladman (2007); Jordan & Bakos (2008) 


general relativistic effect 


< 2 [ms year -1 ] c 


Heyl & Gladman (2007); Jordan & Bakos (2008) 


stellar quadrupole moment 


< 0.004 [ms year -1 ] c 


Heyl & Gladman (2007); Jordan & Bakos (2008) 



a The case having a Jovian-mass secondary planet with an orbital period of ~ 3 years, which corresponds to the actual 
observational period. This produces the RV amplitude of ^21 m s -1 which corresponds to the upper limit placed from the 
RV data. 

b Adopted a value estimated for OGLE-TR-56b, whose properties are similar to the ones of the WASP-5b system, as an 
approximation. 

c Equations of (5), (1), and (3) of Jordan &i Bakos (2008) arc used for calculating precession rates due to tidal deformations, 
the general relativistic effect, and stellar quadrupole moment, respectively, and an equation (23) of Heyl & Gladman (2007) 
into which the rates are substituted is used for calculating the transit timing delays. Eccentricity of the primary planet of 
0.032 is adopted which is 3 ct upper limit derived in this work. 
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Fig. 1. Systematic-corrected light curves of WASP-5b (black points) and their residuals (gray points) 
from the best fit model (solid lines) which was determined by the simultaneous fit for seven transit light 
curves with their normalized error bars taking red noises into account. 0.968 is added to each residual 
light curve for display. E numbers stand for the corresponding transit epochs. 
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Fig. 2. The difference (lower panel) between the light curve corrected using coefficients derived by fitting 
only OOT part (open circles in upper panel) and the one corrected using coefficients derived by fitting 
overall light curve with the transit-model parameters fixing at the values presented in S09 (filled circles 
in upper panel). The initial error bars estimated from the equation (1) arc ovcrplottcd on the differential 
data. 
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Fig. 3. Phase folded, 120-s binned light curves of WASP-5b from the Danish telescope (top) and the B&C 
telescope (second) are shown. Black lines show the best-fit models for the Danish data (results for the 
case using the quadratic limb-darkening law presented in S09) and for the B&C data (results for a joint 
fit to the seven light curves). Third and bottom plots show their residuals (0.965 and 0.970 are added, 
respectively, for display). The rms scatters are 0.00046 and 0.00081, respectively. 
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Fig. 4. Transit Timing residuals of WASP-5b using a new ephemeris determined in section 4. A filled 
square is the Euler data, an open square is the FTS data, two triangles are the Danish data, and seven 
filled circles are the B&C data. 
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Fig. 5. Estimated 3 a upper limits on mass of the hypothetical secondary planet as a function of or- 
bital-period ratio (black thick-solid line), assuming two planets arc co-planer and have initially circular 
orbits. The dotted line represents upper limits on the secondary mass estimated from the RV data, which 
corresponds to a line causing the RV amplitude of 21 m s _1 due to the secondary body. The thick-dashed 
line shows boundary of Hill-stable region calculated from Barnes & Greenberg (2006). Gray vertical-dashed 
lines represent corresponding MMRs. 
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Fig. 6. The best-fit mass regions of the hypothetical secondary planet as a function of orbital period 
ratio (thin-solid lines in upper panels) and their % 2 values (lower panels), for the case of e 2 = (left) and 
0.1 (right). The gray thick lines are same as the black thick line in figure 5. The locations of four example 
TTV models shown in figure 8 are marked as star symbols and alphabets of (a), (b), (c), and (d). 




1.9 1.95 2 2.05 1.9 1.95 2 2.05 

P 2 /P, P 2 /P, 

Fig. 7. Zooms around the 2:1 MMR of figure 6. The pitches of the "observed" period ratio are enlarged 
in the vicinity of the 2:1 MMR, because the strong resonance much excites the orbital eccentricities of the 
two planets resulting in the jumps at a certain period-ratio range. 
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(a) P2/P1 = 0.499475 M2 = 0.003 [Mjup] e2 = 0.0 
X 2 = 4-35 



(b) P2/P1 = 2.062011 M2 = 0.093 [Mjup] e2 = 0.0 
X= 5.01 




(c) P2/P1 = 2.041586 M2 = 0.011 [Mjup] e2 = 0.1 
X 2 = 6.65 



400 
epoch 

(d) P2/P1 = 3.988897 M2 = 0.053 [Mjup] e2 = 0.1 
X 2 =20.57 




400 
epoch 



Fig. 8. Four example TTV models (black curves) with the observed TTV data (legends are same as 
figure 4) (top panels) and their residuals (bottom panels). The explanations for these models appear in 
the text. 



30 



200 



(e) P2/P1 = 3.988897 M2 = 0.053 [Mjup] e2 = 0.1 

X 2 =20.57 



1.2m Euler 
2.0m FTS 
1.5m Danish 
0.6m B&C 



100 - 







100 



■200 
200 






■100 
■200 







200 



400 
epoch 



600 



800 



J 1 

3 


1 

, * - 


1 1 1 1 

f ii i 


i 

I 


-I] 
1 i 


« 


i 1 


► - 



800 



